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Abstract. We study the mathematical model of thermoacoustic tomography in media with a 
variable speed for a fixed time interval [0, T] so that all signals issued from the domain leave it after 
time T. In case of measurements on the whole boundary, we give an explicit solution in terms of a 
Neumann series expansion. We give almost necessary and sufficient conditions for uniqueness and 
stability when the measurements are taken on a part of the boundary. 



1. Introduction 

In thermoacoustic tomography, a short electro-magnetic pulse is sent through a patient's body. 
The tissue reacts and emits an ultrasound wave from any point, that is measured away from 
the body. Then one tries to reconstruct the internal structure of a patient's body form those 
measurements, see e.g, El IHl [HI E2]- For more detail, an extensive list of references, and the 
recent progress in the mathematical understanding of this problem, we refer to [H |H El 13 [121 [H] • 
Both constant and non-constant sound speeds have been studied and naturally, the results are more 
complete in the constant speed case. 

The purpose of this work is to study this problem under the assumption of a variable speed. We 
will actually formulate the problem in anisotropic media. Let 5 be a Riemannian metric in R", 
let o be a vector field, and let c > 0, g > be functions, all smooth and real valued. Assume for 
convenience that g is Euclidean outside a large compact, and c — l = q = a = there (since we 
work with f in a fixed interval, by the finite speed of propagation, this assumption is not essential). 
Let P be the differential operator 

Let u solve the problem 

{d'^ + P)u = in(0,r)xR", 

u\t=o = f, 

dtu\t=o = 0, 

where T > is fixed. 

Assume that / is supported in 0, where fl C R" is some smooth bounded domain. The mea- 
surements are modeled by the operator 

(3) A/ := n|[o,T]x90- 

The problem is to reconstruct the unknown /. 

The presence of the magnetic field {aj} is perhaps of no interest for applications but it does not 
cause any additional difficulties. 
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If T = oo, then one can solve a problem with Cauchy data at t = oo (as a limit), and boundary 
data h = A/. The zero Cauchy data are justified by local energy decay that holds for non-trapping 
geometry, for example (actually, it is always true but much weaker and not uniform in general). 
Then solving the resulting problem backwards recovers /. Now, based on that, one can show 
that for a fixed T, one can still do the same thing with an error e(T) ^ 0, as T ^ cxd. This is 
known as the time reversal method. In the non-trapping case, n odd, the error is uniform and 
e(T) = ©(e""^/*^). There is no good control over C though. Error estimates based on local energy 
decay can be found in [8] , see also Corollary [TJ Other reconstruction methods have been used as 
well, see, e.g., P] for a discussion, and they all use measurements for all t in the variable coefficients 
case, i.e., T = oo; and they are only approximate for T < oo with an error depending on the local 
energy decay rate. Of course, if n is odd and P = —A, any finite T > diam(J7) suffices by the 
Huygens' principle. 

We refer to Section |3] for a discussion of uniqueness results. 

In this paper, we want to study what happens when T < oo is fixed, greater than the length 
of the longest geodesic in (thus the metric c~'^g is assumed to be non-trapping). In case of 
measurements on the whole boundary, our main result is that the problem is Fredholm, uniquely 
solvable, and can be solved explicitly with a Neumann series expansion. In case of partial data, in 
Section [3] we give an almost necessary and sufficient condition for uniqueness, and another almost 
necessary and sufficient condition for stability. In Proposition [3] we characterize A as a sum of two 
Fourier Integral Operators with canonical relations of graph type. 



Notice first that P is formally self-adjoint w.r.t. the measure c ^d Vol, where d Vol(x) = ^det g dx. 
Given a domain U, and a function u{t,x), define the energy 



where Dj = -id/dx^ + aj, D = (L>i, . . . , L>„), \Du\'^ = g'^ {Diu){Dju), and d Vol(x) = (det g)^/^dx. 
In particular, we define the space Hd(U) to be the completion of C^{U) under the Dirichlet norm 



It is easy to see that HoiU) C H^{U), if U is bounded with smooth boundary, therefore, H£){U) 
is topologically equivalent to Hq{U). li U = R", this is true for n > 3 only, if g = 0. By the finite 
speed of propagation, the solution with compactly supported Cauchy data always stays in even 
when n = 2. The energy norm for the Cauchy data (/, ip), that we denote by || • ||-^ is then defined 



2. Complete data 




(4) 




by 




The wave equation then can be written down as the system 



(6) 



= Pu 
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where u = {u,ut) belongs to the energy space 7i. The operator P then extends naturally to a 
skew-self adjoint operator on Ti.. In this paper, we will deal with either U = R" or U = 0,. In the 
latter case, the definition of Hu^U) reflects Dirichlet boundary conditions. 

One method to get an approximate solution of the thermoacoustic problem is the following time 
reversal method, that is actually used in a modified form, see the comments below. Given h, let vq 
solve 



(7) 



(a| + P)vo = in (0, T) X Q, 

^o|[0,T]xaQ = h, 

Vo\t=T = 0, 

dtVo\t=T = 0. 



Then we define the following "approximate inverse" 

AqH := fo(0, •) in 

Then AqA/ is viewed as an approximation to /. As we mentioned above, that is actually true 
asymptotically as T ^ oo, with the modified version of the time reversal method described below, 
(see [8]) but T is fixed in our analysis. 

In this form, the time reversal method has the following downside: h may not vanish on {T} x dQ, 
therefore the mixed problem above has boundary data with a possible jump type of singularity at 
{T} X (the compatibility conditions might be violated). That singularity will propagate back 
to t = and will affect vq, and then vq may not be in the energy space. The operator AqA may 
fail to be Fredholm or even bounded then, and in particular AqA/ might be more singular than /. 
For this reason, h is usually cut off smoothly near t = T, i.e., h is replaced by x(^)^(^)2;), where 
X S C°°(R), X = for t = r, and x = 1 in a neighborhood of (— cjo, T(i7)), see e.g., \Bj Section 2.2]. 

We will modify this approach in a way that would make the problem Fredholm, and will make 
the error operator a contraction. To this end, we proceed as follows. Given h (that eventually will 
be replaced by A/), solve 



{df + P)v = in(0,r)xJ^, 

(8) ' 



v\[o,T]xdn = h-, 
v\t=T = (p, 
dtv\t=T = 0, 



where (p solves the elliptic boundary value problem 

(9) P0 = O, (Plan = h{T,-). 

Since P is a positive operator, is not a Dirichlet eigenvalue of P in $7, and therefore ^ is uniquely 
solvable. Note that the initial data at t = T satisfy compatibility conditions of first order (no jump 
at {r} X (90). Then we define the following pseudo-inverse 

(10) Ah:=v{0,-) inn. 

The operator A maps continuously the closed subspace of H^{[0,T] x (9J7) consisting of functions 
that vanish at t = T (compatibility condition) to if^(r2), see |T3j. It also sends the range of A to 
HQ{n) = H£){n), as the proof below indicates. 

In the next theorem and everywhere below, r(J7) is the supremum of the lengths of all geodesies 
of the metric c~'^g in Cl. Also, dist(x, y) denotes the distance function in that metric. We then call 
{n, c~'^g) non-trapping, if T(J7) < oo. 
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Theorem 1. Let {^D,,c~'^g) be non-trapping, and let T > T(r2). Then AA = Id — K, where K 
is compact in H^^Q), and \\K\\Ho(,i^) ^ ^- -^^ particular, Id — K is invertible on Hu^O,), and the 
inverse thermoacoustic problem has an explicit solution of the form 



(11) 

Proof. Let w solve 
(12) 



f=f2 K'^Ah, h := A/. 



m=0 



w\[o,T]xdn 

w\t=T 
Wt\t=T 





0, 

u\t=T - 
Ut\t=T, 



in (0,r) X 



where u solves ([2]) with a given / S Hd. Let v be the solution of ([8| with h = A/. Then v + w 
solves the same initial boundary value problem in [0, T] x J7 that u does (with initial conditions at 
t = T), therefore u = v + w. Restrict this to f = to get 

f = AAf + w{0,-)- 

Therefore, 

Kf = w{0,-)- 

In what follows, (•, ■)hd{^) inner product in H£){Q,), see (j4]), applied to functions that belong 

to II^{i}) but maybe not to H£){^) (because they may not vanish on d^). Set u'^ := u{T,-). By 
^ and the fact that = (/) on d^, we get 



(n^ - cP, 



Then 



\Hnm 



\u 



T\\2 



< u 



T\\2 



2 



Therefore, the energy of the initial conditions in (12) satisfies the inequality 

(13) En{w,T) = \\u^ - m^^^) + ll^f Ili2(f,) < En{u,T). 
Since the Dirichlet boundary condition is energy preserving, we get that 

En{w,{)) = En{w,T) < En{u,T) < Enn{u,T) = En{u,0) = 
In particular, 

(14) \\Kf\\l^^^^<En{w,0)<\\ffH^^^y 
We show next that actually the inequality above is strict, i.e., 

(15) \\Kf\\H,in)<\\f\\HMn), / / 0. 

Assume the opposite. Then for some / 7^ 0, all inequalities leading to (14) are equalities. In 
particular, Eq{w,T) = E-^n(u,T). Then 

n(T, x) = 0, for x il. 

By the finite domain of dependence then 

(16) u{t,x)=0 when dist(a;,J7) > |r - t|. 
One the other hand, we also have 

(17) u{t,x) = when dist(x,0) > \t\. 



THERMOACOUSTIC TOMOGRAPHY 



5 



Therefore, 

(18) u{t, x) =0 when dist(a;, dn) > T/2, -T/2 < t < 3T/2. 

Since u extends to an even function of t that is still a solution of the wave equation, we get that 
(18) actually holds for \t\ < 3T/2. Then one concludes by Tataru's theorem, see Theorem |4j that 



It = on [0,T] X Q, therefore, / = 0. We refer to [5| for a similar argument. Note that the time 
interval here is actually larger than what we need for the uniqueness argument, see also Theorem [2] 
and Corollary [2] below. 

We will show now that K is compact. Since T > T(J7), all singularities starting from 0, leave 
ai t = T. Therefore, u(T,-) and ut(T,-), restricted to Cl, are C°° . Moreover, considered as linear 
operators of /, they are operators (FIOs, actually) with smooth Schwartz kernels. Then so is (p, 
see by elliptic regularity. Therefore, the map Hd{Q,) 9 / i— > u{T, •) — G H£,{^) is compact 
because it is an operator with smooth kernel on 0. Next, the map Hd(0,) 9 / i— > ut(T, •) G Ho{^) is 
compact as well. Since the solution operator of ( 12 ) from t = T to t = is unitary in i/£i(J7)©L^(J7), 



we get that the map Hd{^) B f ^ w(0, •) G Hd^Q) is compact, too, as a composition of a compact 
and a bounded one. 
Now, one has 



(19) \\Kf\\H,in)<V>^i\\f\\Hnm, //O, 



where Ai is the largest eigenvalue of K*K. Then Ai < 1 by (15). □ 



Remark 1. Although we proved that K is compact, we did not show that K is smoothing of 1 
degree. Actually, we showed that is a composition of a smoothing operator and a bounded one. 



To make K smoothing, we need to modify the initial condition for wt(T, •) in (12), as we did it for 
w{T, •), so that it would satisfy the compatibility at {T} x 90 (no jump there, i.e, wt{T, •) G Hq{^1)). 
That will put {wiT, ■),wt(T, •)) in the domain of the generator of the solution group, in other words, 
{wt, Pw(T, •)) would be in the energy space. Then the same would be true for Pw{0, •) = —PKf, 
hence Kf G //^(il). Then we get a Fredholm problem again but the norm of K may not be less 
than 1 (that still might be true in a suitable norm). In any case. Id — K will be invertible. One can 



also modify the initial data at t = T in (12) to satisfy even higher order compatibility condition. 



and that will increase the smoothing properties of K. 

Remark 2 . The smoothness requirements on the coefficients of P can be relaxed to require smooth- 
ness of a finite degree. All we need, besides a well posed problem in the energy space, is a propaga- 
tion of singularities result with a gain of smoothness on t = T enough to guarantee compactness of 
K; and Tataru's uniqueness theorem in that case. We will not pursue this for the sake of simplicity 
of the exposition. 

The proof of Theorem [T] provides an estimate of the error in the reconstruction if we use the first 



term in (11 ) only that is Ah. It is in the spirit of and relates the error to the local energy decay. 



as can be expected. 
Corollary 1. 

11/ - AA/y^^n) < (fl^) ' WfWHnin), V/ G H^^n), / / 0, 
where u is the solution of 

Note that the / — AAf = Kf, and the corollary actually provides an upper bound for ||-fC/||. 
The estimate above also can be used to estimate the rate of convergence of the Neumann series (llll) 
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when we have a good control over the uniform local energy decay from time t = to time t = T. The 
estimate holds even without the non-trapping condition and for any T > but Efi{u,T)/E^{u,0), 
that is always less or equal to 1 , can be guaranteed to have a uniform upper bound less than 1 for all 
/ only when T > T(0,); then the operator norm of K is less than 1, as well. If T{Q) /2 < T < T(0), 
we can only say that ||-ftr/|| < ||/|| for any /, see Corollary [2] below but that does not necessarily 
imply that ||i^|| < 1. If T < T{Q)/2, then there is always / so that that quantity equals 1 by a 
trivial domain of dependence argument. 

3. Incomplete data 

The case of partial measurements has been discussed in the literature as well, see e.g..[T2 l l23 l [2i]. 
One of the motivations is that in breast imaging, for example, measurements are possible only on 
part of the boundary. For simplicity, we assume in this section that P = —A outside in particular 
c = 1 and g is Euclidean outside 

(20) c{x) = 1, gij{x) = 5ij, for x ^ Q. 

All geodesies below are related to the metric c~^g. 
Let r C (90 be a relatively open subset of dU. Set 

(21) g := {{t,x); X eT, <t < s{x)} , 

where s is a fixed continuous function on T. This corresponds to measurements taken at each x gT 
for the time interval < t < s{x). The special case studied so far is s{x) = T, for some T > 0; 
then g = [0, T] x T. 

We assume now that the observations are made on g only, i.e., we assume we are given 

(22) Af\g, 

where, with some abuse of notation, we denote by A the operator in ([s]), with T = oo (that 
actually can be replaced by any upper bound of the function s). Then we want to know under 
what conditions one can recover /, and when that recovery is stable. 

Uniqueness and reconstruction results in the constant coefficients case based on spherical means 
were known for a while, see e.g., the review paper [12\. If P = —c'^{x)A, and g = [0, T] x 50, Finch 
and Rakesh [4^ have proved that A/ recovers / uniquely as long as T > T(0). A uniqueness result 
when r is a part of 50 in the constant coefficients case is given in |3], and we follow the ideas of 
that proof below. The Holmgren's uniqueness theorem for constant coefficients and its analogue 
for variable ones, see Theorem |4] below, play a central role in the proofs that suggests possible 
instability without further assumptions, see also the remark following Theorem |3] below. Stability 
of the reconstruction when P = — A and T = oo follows from the known reconstruction formulas, 
see e.g., [12]. In the variable coefficients case, stability estimates as T — > oo based on local energy 
decay have been established recently in [8j. When T is fixed, there is the general feeling that if one 
can recover "stably" all singularities, and if there is uniqueness, there must be stability (although 
this has been viewed from the point of view of integral geometry, see also Section |4| . We prove this 
to be the case in Theorem [sj and we use analysis in fT6] , as well. 

We present some heuristic arguments for our main assumption below. We will restrict / below 
to a class of functions with support in some fixed compact /C C O. Intuitively, to be able to recover 
all / supported in /C, we want for any x G /C, at least one signal from x to reach g, i.e., we want to 
have a signal that reaches some z £T for t < s{z). In other words, we should at least require that 



(23) 



yx £ IC,3z £ T so that dist(x, z) < s{z), 
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(one may want to have a non-strict equality above but we will not pursue this). In Theorem |4] 
below, we show that this necessary condition, up to replacing the < sign by the < one, is sufficient, 
as well. 

If we want that recovery to be stable, we need to be able to recover all singularities of / "in a 
stable way." By the zero initial velocity condition, each singularity (x, ^) splits into two parts, see 
Proposition [3] below: one that starts propagating in the direction ^; and another one propagates in 
the direction — ^. Moreover, neither one of those singularities vanishes at t = (and therefore never 
vanishes), they actually start with equal amplitudes. For a stable recovery, we need to be able to 
detect at least one of them, in the spirit of |16] . i.e., at least one of them should reach Q. This in 
particular allows us to reduce T by half in the full boundary data case, i.e., when Q = (0, T) x dQ, 
one can choose 

(24) T>T{n)/2, 

and still hope that a stable recovery is possible. In the general case, define t±{x, ^) by the condition 

T±ix,C) = max (r > 0; -fx,i{±r) £ 0,) . 

Based on the arguments above, for stable recovery we should assume that Q satisfies the following 
condition 

(25) V(x,^) G 5*/C, (tct(x, ^), 7^^^(ro-(x, ^)) G Q for either o" = + or cr = — (or both). 



Compared to condition (23), this means that for each x G /C and each unit direction ^, at least one 
of the signals from (x, ^) and (x, — ^) reaches Q. This condition becomes necessary, if we replace Q 
by its closure above, see Remark |4} In Theorem |3] below, we show that it is also sufficient. 

3.1. Uniqueness. We have the following uniqueness result, that in particular generalizes the result 
in [3] to the variable coefficients case. 

Remark 3. Note that we do not need the geodesic flow to be non-trapping in this theorem since 
(23) is a condition on a subset of the geodesies only. 

Theorem 2. Let P = —A outside 0,, and let dO, be strictly convex. Then under the assumption 
[2^, ifAf = Oong for f G Hd{^) with supp/ C IC, then / = 0. 

Proof. We follow the proof in [3], where g is Euclidean everywhere, and T = 00 (actually, it is easy 
to see there that T can be any number larger than T(Q)). We preserve the notation of |3] as much 
as possible. 

Recall that dist(x, y) is the distance in the metric c~^g. Let d(x, y) be the (Euclidean) distance 
in R" \ 17 defined as the infimum of the Euclidean length of all smooth curves in R" \ joining 
X and y. The function d is Lipschitz continuous, see |3j. Let Er{x) be the ball with center x and 
radius r > in that metric. Then in [3 Proposition 2] , Finch et al. proved the following domain of 
dependence results for solutions vanishing on a part of dQ. 

Proposition 1 ([3j). Let 0, be an open bounded connected subset 0/ R" with a smooth boundary. 
Suppose u is a smooth solution of the exterior problem 

utt- Au = 0, t G R; X G R" \ n, 
u = h on R X dU,. 

Choose p ^ 0,, and to < ti. If u(tQ, •) = ut{tQ,-) =0 on Etj^-toip)j ^ ^-^ zero on 

{{t, x); X G d^l, to < t < ti, d{x,p) < ti — t} , 
then u(t, p) = (t, p) = for all t £ [to, ^i] • 
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Let A/ = on ^, with / as in the theorem, and let u be the corresponding solution of 
Fix a point xq G fC. We will show that / = near xq. By (23), there is p G Oil so that 
dist(xo,p) < s(p); then {s{p),p) gQ. Let < p <C 1 be such that [0, s{p) - p] x {Ep{p) n dVL) C Q, 
and dist(xo, q) < s{q) — p, \/q £ Ep{p) fl dU. We can therefore assume that 

(26) g = [0, T] X r, where T = Ep{p) n 90, 
and 

(27) dist(xo,g)<T G L. 
The first step of the proof if to show that 

(28) / = inBp{p), 



where Bp{p) is the ball in the metric g with center p and radius p. The proof of (28) is the same as 
in [3] with taking extra care about the range of the t variable. Indeed, notice first that u solves the 
wave equation in R" \ Q with zero Cauchy data there, and boundary data h = ti|R, X(9Q vanishing 



on Q, see (26). Fix a small neighborhood U oi p outside By ( |27[ ) and the finite domain of 
dependence result in Proposition [T| we get n = on {—p + e, p — e) x U, where < e ^ 0, when 
the size of U tends to 0. 

Next, u solves the wave equation in the whole space, and can be extended (as a solution) as an 
even function of t. Therefore, by the unique continuation principle, see Theorem |4j we get (28). 

The next step is to iterate this argument and to prove that / = near xq. This would follow 
from the following property that we prove next: For some o" > independent of p, we have 

(29) / = mBr{p),r>p =^ f = in S^i„|^+,,r}(p). 



The reason we did not just replace the minimum above with p+a is that we apply (29 ) consecutively 
several times; at each step we gain a, and we would like to make the radius equal to T. The last 
step needed for that might be smaller than a though, and (26), (27) pose a restriction on how far 
we can go. 

Relation (29) follows from the following. 

Lemma 1. Assume that supp/ C K = 0, \ Br{p) with some r > p. Let 5 = (i\st{E p{p) ^ K) . Then 
f = in B^i^{p+s,T}{p)- 

We prove Lemma [l] below. Let a be the supremum of the distance dist(p, q), g G F. Since dO, is 
strictly convex, a < p. Indeed, a is actually the maximum of those distances, if we replace F by 
the compact F. Then a = dist(p, qo) for some qo £ T. Because of the strict convexity the latter is 
the length of the shortest geodesic on connecting p and go- If we assume that a = p, then that 
geodesic will be a minimizing curve for c~'^g as well, therefore it will be a geodesic for that metric. 
That is impossible because for p <^ 1, there is unique minimizing geodesic connecting p and go, and 
that geodesic cannot be on dfl. 

The following lemma generalizes |3l Propositon 5] to the current setting. We refer to Fig. 1 that 
is similar to Fig. 2.5 there for better understanding of the lemma and its proof. 

Lemma 2. dist{K, Ep{p)) is the length of some geodesic segment joining a point A G K and a 
point B gT. 

The proof is provided below, and we continue with the proof of Theorem [2] By Lemma [2] 5 is 
the length of the geodesic segment connecting A and B, as in the lemma. Then 

p + 5 = p + dist{A,B) =dist{A,B) +dist{B,p) + {p - dist{B,p)) 

^^^^ > dist{p,A) + {p-dist{B,p)) >r + {p-a). 
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Figure 1. Illustrates Lemma [T| One can also show that A G dil H dBr{p). 



Note that a := p — a > is independent of r. This proves the property (29), and therefore, the 
theorem. □ 

It remains to prove the two lemmas above. 

Proof of Lemma We will provide a proof that is different and shorter than that in [3] . Since 
dist {K, Ep{p)) is the distance between two compact sets, there is A G K and B G Ep{p) so that 
dist{K, Ep{p)) = dist(^, i?). By the Hopf-Rynow theorem, there is a geodesic 7 connecting A, B 
so that Aist{A,B) is the length of 7. Clearly, B belongs to dEp{p) that consists of two parts: the 



first one that we denote by dEp^^{p), that is outside and the second one is T, see (26) We will 



show first that B must belong to the second one. Assume the opposite. Then 7 intersects 90 once 
(because of the strict convexity) at some point C T because if C S T, then we would have C = B. 
Then the segment Ci? of 7 is a straight line segment outside Ep{p), see Figure 2. 




Figure 2. 

Let c be the minimizing curve for the metric d, lying outside O, that connects B and p. It is easy 
to see (see [3]) that c exists and consists of a straight line segment ci = BD between C and some 
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D G dEp{p), and a geodesic C2 on d^, possibly reduced to a point, so that ci and C2 are tangent to 
each other and to at their common point that we denote by D. Note that dEp^^{p) is an open 
surface, therefore B ^ D. Then the curve CB U BD locally minimizes the lengths of all curves 
connecting C and D with the property that they consist of a curve outside Ep{p) Ui7 connecting C 
and some B' G dE^^^{p) close to B] and then another curve, outside O but inside Ep(p), connecting 
B' to D. Then CB U BD must be a straight line segment; otherwise we can make it shorter by 
an arbitrary small perturbation, and that would contradict the minimizing property above. That 
segment is tangent to dQ. By the strict convexity of J7, it cannot have two common points C and 
D with do,. This contradiction shows that B £ T, and this proves the second statement of the 
lemma. □ 

Proof of Lemma [1} Roughly speaking, the idea of the proof is that we can apply the arguments at 
the beginning of the proof of the theorem by shifting the initial moment form t = to t = 5. 
First, by the definition of 5 and the standard domain of dependence argument, 

(31) u = ut = on [-6, 6] x Ep{p). 

Let U he a small enough neighborhood of p in Ep{p). li 5 + p <T, by the domain of dependence 



argument for the exterior problem [3l Proposition 2], in view of (26), (27), u = on [5,6 + p 
o(l)] X [/, where by o(l) we denote terms tending to when the size of U tends to 0. If 5 + /) > T, 
then we can prove that only in the time interval [5, T — o(l)]. Therefore, in both cases, the time 
interval is [5, min{5 + p, T} — o(l)]. Since u extends as an even solution in the t variable, we get 
that ti = for \t\ < min{5 + p,T} — o(l), x £ U. Then from the unique continuation result in 
Theorem|4j u\t=o, and therefore / vanishes in -Bmin{<5+p,T}-o(i)(p)- Letting U tend to p, we get that 
f = in B^i^{S+p,T}{p)- □ 

It is probably worth mentioning that we actually proved the following result about partial re- 
covery given insufficient information. 

Proposition 2. Let P = —A outside J7, and let dQ be strictly convex. Assume that Af = on Q 



for some f G HD{i}) with supp/ C with Q as in (20) that may not satisfy (23). Then f = in 
W , where 

W := {x £0,; G r so that dist(x, z) < s{z)] . 
Moreover, no information about f on 0,\W is contained in Af\g. 

3.2. Stability. In this section, we use tools from microlocal analysis. We refer, for example, to 
|20j for an introduction to the theory of pseudo-differential operators (\I'DOs) and to \2T\ |2] for the 
theory of Fourier Integral Operators (FIOs). 

We now consider the situation where Af is given on a set Q satisfying (p5|). Since /C is compact 



and Q is closed, one can always choose Q' (s^ Q that still satisfies (25). Fix x £ C'o°([0)^] ^ 90,) so 
that suppx C G and x = 1 on G'. The measurements are then modeled by X-^fi which depends on 
Af on Q only. 

We start with a description of the operator A that is of independent interest as well. In the next 
proposition, we formally choose T = oo. 

Proposition 3. A = A_|_-|-A_, where A± : Cq°{0) —>■ C°°((0, oo) x dO) are elliptic Fourier Integral 
Operators of zeroth order with canonical relations given by the graphs of the maps 

(32) {y,0 ^ ir±{y,0,iy,±d^±{y,0),\^\,iy,±i{r±{y,0)) , 

where \(,\ is the norm in the metric c~'^g, and the prime in j' stands for the tangential projection 
of 7 on TdO. 



THERMOACOUSTIC TOMOGRAPHY 



11 



Proof. This statement is well known and follows directly from [2 , for example. We will give more 
details that are needed just for the proof of this proposition in order to be able to compute the 
principal symbol in Theorem [3j 

We start with a standard geometric optics construction. Fix xq G O. In a neighborhood of 
(0,xo), the solution to ([8| is given by 

(33) uit, x) = (27r)-" [ e''^'^(*'"'«)a.(x, ^, t)f{C) d?, 

modulo smooth terms, where the phase functions (j)± are positively homogeneous of order 1 in ^ 
and solve the eikonal equations 

(34) =F 94(/>± = |d^(/)±|, (t)±\t=Q = X ■ 

while a± are classical amplitudes of order solving the corresponding transport equations, see [21 
p. 128] or |2H eqn. (VI. 1.50)]. In particular, a± satisfy 

a+ + a_ = 1 for t = 0. 

Since dt4'± = =F'^ for t = 0, and = for t = 0, we also see that 

a+ = a_ for t = 0. 

Therefore, a_|_ = a_ = 1/2 at t = 0. Note that if P = A, then (j)± = x ■ £,^t\^\, and 0+ = a_ = 1/2. 
The principal term a"^ of a± ~ X]j>o '^L '''' satisfies the homogeneous transport equation 

(35) [dt - c^g'Hd,,ct)±)d,, + C±) a± = 0, a±\t=o = 1/2, 

where Cj depend on the coefficients of P and on (p±, see |21[ eqn. (VI. 1.49)]. 

By the stationary phase method, singularities starting from (x, ^) G WF(/) propagate along 
geodesies in the phase space issued from (x, ^), for a = +. i.e., they stay on the curve {'yx,({t),'yx,({t))', 
and from (x, — ^), for a = —, i.e., they stay on the curve (jx-^it) , ix,-^{'t)) . This is consistent with 
the general propagation of singularities theory for the wave equation because the principal symbol 
of the wave operator — c^lCls has two roots r = ibc|^|g. 

The construction is valid as long as the eikonal equations are solvable, i.e., along geodesies 
issued from (x, it^) that do not have conjugate points. Assume that WF(/) is supported in a 
small neighborhood of (xo,Co) with some 0- Assume first that the geodesic from (xo,Co) with 
endpoint on 917 has no conjugate points. We will study the cr = + term in ( |33| ) first. Let (/)b, ab 
be the restrictions of (/>+, a^, respectively, on R x dO,. Then, modulo smooth terms, 

(36) A+/ := u+{t,x)\nxan = (2^)^" / ei*'^(*'^'«)ab(x, t)/(0 d^, 



where u+ is the a = + term in (33). Set to = t+(xo,Co), yo = 7xo,?o(*o), % = 7a;o,5o(*o); in other 



words, (yO)%) is the exit point and direction of the geodesic issued from (xo,^o) when it reaches 



dQ. Let X = (x',x"') be boundary normal coordinates near yo- Writing / in (36) as an integral, we 
see that ( pG) ) is an oscillating integral with phase function <I> = (/>+(t, x', 0,5,) ~ V ' ^- Then (see [21] . 
for example), the set T, := = 0} is given by the equation 

It is well known, see e.g.. Example 2.1 in [211 VI. 2], that this equation implies that (x',0) is the 
endpoint of the geodesic issued from (y,^) until it reaches the boundary, and t = T^{y,^), i.e., t 
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is the time it takes to reach dO,. In particular, S is a manifold of dimension 2n, parametrized by 
(y,^). Next, the map 

(37) E 3 {y, t, x', ^ ' — ' {v^ 2;', dt(l)+, 

is smooth of rank 2n at any point. This shows that $ is a non-degenerate phase, see [21 \ VIII. 1], 
and that / 1-^ A^f is an FIO associated with the Lagrangian given by the r.h.s. of ( |37[ ). The 
canonical relation is then given by 



Then ( 32 ) follows from the way is constructed by the Hamilton- Jacobi theory. The proof in the 
cr = — case is the same. 

The proof above was done under the assumption that there are no conjugate points on 7j;o,^o(^)' 
< t < r+(yO)^o)- To prove the theorem in the general case, let ti G (0, r+(yo, Co)) be such that 
there are no conjugate points on that geodesic for ti < t < T+(yO)Co)- Then each of the terms in 
( |33[ ) extends to a global elliptic FIO mapping initial data at t = to a solution at t = ti, see e.g., 
|2]. Its canonical relation is the graph of the geodesic flow between those two moments of time 
(for a = +, and with obvious sign changes when a = —). We can compose this with the local FIO 
constructed above, and the result is a well defined elliptic FIO of order with canonical relation 
321). □ 



Choose and fix T > supp s, see (21). Let A be the "back-projection" operator defined in and 



(10). Note that A is always applied to below, therefore (/> = in this case. 



Theorem 3. Ax-A is a zero order classical ^DO in some neighborhood of fC with principal symbol 

\ {X{lx,d^+{x, m + X(7x,c(r-(x, e)))) . 



If Q satisfies (25), then 

(a) AxA is elliptic, 

(b) AxA is a Fredholm operator on H£,{1C), and 

(c) there exists a constant C > so that 

(38) ll/IU,(yc) < CWKfWmigy 



Remark 4. By |16| Proposition 3], condition (25), with Q replaced by its closure, is a necessary 
condition for stability in any pair of Sobolev spaces. In particular, c~'^g has to be non-trapping 
for stability. Indeed, then the proof below shows that Ax^ will be a smoothing operator on some 
non-empty open conic subset of r*/C \ 0. 

Remark 5. Note that A : Hd^IC) — > H^{[0,T] x 917) is bounded. This follows for example from 
Proposition |3| 

Proof. We will use the geometric optics construction in the proof of Proposition [3] using the notation 
there. 

To construct a parametrix for AxAf, we apply a geometric optic construction again, using the 
two characteristic roots ibc|^|g. The boundary data A_|_/ have a wave front set in a small conic 
neighborhood of ((to, yo)' (1' ^o))- Note that ^ because geodesies issued from /C cannot be 
tangent to 50. Then for the solution w of ([s]) with h = A4./, we can apply the geometric optics con- 
struction above, but now with initial condition on R x dQ, to get two types of singularities starting 
from that one. The first one propagates along the geodesies close to 7^0 ,5o ™ opposite direction. 
The second one propagates along the geodesic close to the one issued from ((to, yo), (^', —V^))-, that 
is transversal to 00. This ray is in fact a reflected 7xo,^o- the propagation of singularities 
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results, those singularities stay on those geodesies until they reach dO, again, then reflect by law of 
geometric optics, etc., i.e., they propagate along the broken geodesies issued from a neighborhood 
of that point. Near t = T however, the solution to (|8|, where = 0, is zero because we have zero 
Cauchy data, and h = X-^f = for t close to T. This shows that the second types of singularities 
do not exist; and in our parametrix construction, we need to work with the first one only. 

We look for a parametrix of the solution of the wave equation ([8| with zero Cauchy data at t = T 
and boundary data X-^+f the form 

v{t,x) = {27Ty^ [ e*'-'«)6(x,e,t)/(e)de. 



The arguments above show that (p = <j)^. Next, for x G dO,, we have b = x^- We need to find h at 
t = 0. The amplitude h satisfies the same transport equation as in the proof of Proposition [3] but 
with initial condition at R x 90. In particular, it is a classical amplitude of order 0. Let 6o be its 



principal part. Then 6o satisfies (|35|), also satisfied by that is a linear homogeneous ODE along 



the bicharacteristic close to (7x0,50) 7x0,^0 )• Therefore, is a linear function of its initial condition 
at R X If we assume for a moment that x = li then we would get 6o = cl'^ , therefore, 6o = 1/2 
for t = 0. Therefore, we get that hQ{x,^)\t=o is given by the value of x/2 at the exit point of 7^;^^ 
on 90 because that value is the initial condition of the transport equation on that bicharacteristic. 

The arguments above reveal the geometry of the singularities but some of them are not needed 
for the formal proof. One can define v as above, localized near the bicharacteristic issued from 
{xq: Co)i and let m+ be the solution of ([s]) with = and h = x^+f ■ Then one easily checks that 
w := u^—v solves the wave equation modulo smooth terms, with smooth boundary condition, and 
that w = near t = T; and is therefore smooth. 

In the same way one treats the a = — term. This proves the theorem assuming no conjugate 
points in $7. 

In the general case, we can apply those arguments step by step, in intervals [0,ti], then [ti,t2]) 
etc., short enough so that there are no conjugate points on the corresponding geodesic segments. 
After the first step, we get {u,ut) at t = ti. Then we construct a parametrix from t = ti to t = t2 
using a new phase function. Note that now, when a = +, for example, ut\t=ti does not vanish 
anymore. On the other hand, {u,ut)\t=ti is Cauchy data of a solution which singularities do not 
travel in two opposite directions, and we will still get one term only, that is an analogue of the 



o" = + one in (33). Then we reach the boundary and apply the result above. Next, step by step, we 
go back to the hyperplane t = 0. An alternative way is to apply the Egorov's theorem from t = 
to t = t, instead of the partition of the time interval, where t is such that there are no conjugate 
points on the bicharacteristic issued from (xq, ■^0) from t to t^{xo, ^0); and on that segment, we use 
the arguments above. 

This proves the first statement of the theorem. 



Parts (a), (b) follows immediately from the ellipticity of AxA that is guaranteed by (25). 
To prove part (c), note first that the ellipticity of AxA and the mapping property of A, see |13] . 
imply the estimate 

\\f\\HniK)<C{\\xAf\\H^ + \\fh2). 

By Theorem [2| and (25), x-^ is injective on H£){}C). By [19, Proposition V.3.1], one gets estimate 
(38) with a constant C > possibly different than the one above. □ 



Corollary 2. Let g be Euclidean outside 0,, and let dO, be strictly convex. Then if Af = on 
[0,T] X dn for some f £ Hd{^), with T > T{Q)/2, then / = 0. 
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4. THERMOACOUSTIC tomography and INTEGRAL GEOMETRY 

If P = —A, and if n is odd, the solution of the wave equation can be expressed in terms of 
spherical means, as it is well known. Then the problem can be formulated as an integral geometry 
problem — recover / from integrals over spheres centered at dfl, with radii in [0, T], and this point 
of view has been exploited a lot in the literature. One may attempt to apply the same approach 
in the variable coefficients case; then one has to integrate over geodesic spheres. This has two 
drawbacks. First, those integrals represent the leading order terms of the solution operator only, 
not the whole solution. That would still be enough for constructing a parametrix however but 
not the Neumann series solution in Theorem [T| The second problem is that the geodesic spheres 
become degenerate in presence of caustics. The wave equation viewpoint that we use in this paper 
is not sensitive to caustics. We still have to require that the metric be non-trapping in some of our 
theorems. By the remark following Theorem |3] however, this is a necessary condition for stability. 



On the other hand, it is not needed for the uniqueness result as long as (23) is satisfied. 
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Appendix A. Unique gontinuation for the wave equation 

We recall here a Holmgren's type of theorem for the wave equation {d^ + P)u = due mainly 
to Tataru. While this theorem is well known and used, and follows directly from the results cited 
below, we did not find it clearly formulated in the literature. 

Theorem 4. Let P be the differential operator in R" as in the Introduction. Assume that u G H^^^ 
satisfies 

{df + P)u = 

in a neighborhood of [—T,T] x {xq}, with some T > 0, xq £ R". Then 

u{t,x) = for |t| + dist(xo, x) < T. 

Proof. If P has analytic coefficients, this is Holmgren's theorem. In the non-analytic coefficients 
case, a version of this theorem was proved by Robbiano [15] with p replaced by Kp with an 
unspecified constant K > 0. It is derived there from a local unique continuation theorem across 
a surface that is "not too close to being characteristic". In [7 , Hormander showed that one can 
choose K = Y^27/23, in both the local theorem [Tj Thm 1] and the global theorem [7, Corollary 7] . 
Moreover, he showed that K in the global one can be chosen to be the same as the K in the local 
one. Finally, Tataru jl7| IT8] proved a unique continuation result that implies unique continuation 
across any non-characteristic surface. This shows that actually K = 1 in Hormander's work, and 
the theorem above then follows from [3 Corollary 7]. □ 
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